/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: RGrid.cxx,v $

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include <iostream>
#include <string.h>
#include <time.h> // for seeding rand()
#include <limits> //for clearing cin stream on bad input data

#include "vtk_include.h"

using namespace std;

vtkDoubleArray *scalars[2];
vtkImageData *idata;
vtkThresholdPoints *tpoints;
vtkPolyDataMapper *idataMapper;
vtkActor *idataActor;
vtkRenderer *renderer;
vtkRenderWindow *renWin;
vtkRenderWindowInteractor *iren;
vtkBMPReader *reader;
vtkLookupTable *lut;
vtkImageRGBToHSV *hsv;

int ****universe; // 4-dimensional pointer -> [2][dimX][dimY][dimZ]
// two 3-diminsional grids, one is read from, other written to
// stores the age of each point (0 is dead)
int n = 1; // read universe
int n2 = 0; // write universe
// 1 or 0, determines which 'universe' is currently being read/written
// also "scalars" borrows this variable for inheritance coloring (same principle)

int bornLow; // min neighbors for dead point to turn on
int bornHigh; // max neighbors for dead point to turn on
int liveLow; // min neighbors for living point to stay on
int liveHigh; // min neighbors for living point to turn on
int maxAge; // maximum age before point dies
int ageOrInherit; // coloring mode: 0 = age, 1 = inheritance
int dimX = 50;
int dimY = 50;
int dimZ = 50;
int timer = 150; // milliseconds between each time step

// used for allocation and indexing into scalars
int totalSize; // total number of points
int layerSize; // total number of points in an Y*Z layer

int useImage; // 1 if input is image, 0 is randomized

// randomize the universe
void randomize()
{
	int i, j, k, r, gridOffset = 0;

	for(i = 0; i < dimX; i++)
	{
		for(j = 0; j < dimY; j++)
		{
			for(k = 0; k < dimZ; k++)
			{
				r = rand();
				if(r%13 == 0) // about half of all 26 neighbors will be alive
				{
					r = rand();
					if(maxAge)
					{
						// set the age between 1 and maxAge
						universe[0][i][j][k] = r%maxAge + 1;
						universe[1][i][j][k] = r%maxAge + 1;
					}

					if(ageOrInherit) // with inheritance coloring, just randomize colors
					{
						// set the scalar so (0 <= scalar < 1) because
						// dead points (scalars = 1) are filtered out by ThresholdPoints
						scalars[0]->SetValue(gridOffset, r/double(RAND_MAX + 1));
						scalars[1]->SetValue(gridOffset, r/double(RAND_MAX + 1));
					}
					else // with age coloring, match the color to the random age above
					{
						scalars[0]->SetValue(gridOffset, double(universe[0][i][j][k] - 1)/maxAge);
						scalars[1]->SetValue(gridOffset, double(universe[0][i][j][k] - 1)/maxAge);
					}
				}
				else
				{
					// set dead
					universe[0][i][j][k] = 0;
					universe[1][i][j][k] = 0;
					scalars[0]->SetValue(gridOffset, 1);
					scalars[1]->SetValue(gridOffset, 1);
				}
				gridOffset++;
			}
		}
	}
}

void inputImage()
{
	double hue, saturation, value;
	int i, j, k, r, gridOffset = 0;

	// use full hue range (0 to 1) when inheritance coloring enabled
	// otherwise the default (blue to red) lookup table is used
	if(ageOrInherit)
	{
		lut = vtkLookupTable::New();
		lut->SetTableRange(0, 1);
		lut->SetHueRange(0, 1);
		lut->SetSaturationRange(1, 1);
		lut->SetValueRange(1, 1);
		lut->Build();
		scalars[0]->SetLookupTable(lut);
		scalars[1]->SetLookupTable(lut);
	}

	for(i = 0; i < dimX; i++)
	{
		for(j = 0; j < dimY; j++)
		{
			hue = hsv->GetOutput()->GetScalarComponentAsDouble(i, j, 0, 0);
			saturation = hsv->GetOutput()->GetScalarComponentAsDouble(i, j, 0, 1);
			value = hsv->GetOutput()->GetScalarComponentAsDouble(i, j, 0, 2);

			//cout << "\n" << i << "," << j << ":" << hue << " " << saturation << " " << value;
			for(k = 0; k < dimZ; k++)
			{
				r = rand()%255;

				// higher value (intensity) increases chance of live points along z axis
				if(r < value)
				{
					if(ageOrInherit)
					{
						// with inheritance coloring, use hue to set scalar
						universe[0][i][j][k] = 1;
						universe[1][i][j][k] = 1;
						scalars[0]->SetValue(gridOffset, hue/255 - 0.001);
						scalars[1]->SetValue(gridOffset, hue/255 - 0.001);
					}
					else
					{
						// with age coloring, map the saturation to nearest age increment, then set color
						universe[0][i][j][k] = saturation*maxAge/255.001 + 1;
						universe[1][i][j][k] = saturation*maxAge/255.001 + 1;
						scalars[0]->SetValue(gridOffset, double(universe[0][i][j][k] - 1)/maxAge);
						scalars[1]->SetValue(gridOffset, double(universe[0][i][j][k] - 1)/maxAge);
					}
				}
				else
				{
					universe[0][i][j][k] = 0;
					universe[1][i][j][k] = 0;
					scalars[0]->SetValue(gridOffset, 1);
					scalars[1]->SetValue(gridOffset, 1);
				}
				gridOffset++;
			}
		}
	}
	
}

int getIntInRange(const char *prompt, int min, int max)
{
	int input;

	while(true)
	{
		cout << prompt;
		cin >> input;

		// bad input
		if(!cin.good())
		{
			cin.clear();
			cin.ignore(numeric_limits<streamsize>::max(),'\n');
			cout << "The value must be an integer.\n";
		}
		else if(input < min)
			cout << "The value must be greater than or equal to " << min << ".\n";
		else if(input > max)
			cout << "The value must be less than or equal to " << max << ".\n";
		else
			return input;
	}
}

void getRules()
{
	reader = vtkBMPReader::New();
	int dim[3];
	char imagePath[128];

	bornLow = getIntInRange("Enter minimum neighbors to be born.\n", 1, 26);
	bornHigh = getIntInRange("Enter maximum neighbors to be born.\n", bornLow, 26);
	liveLow = getIntInRange("Enter minimum neighbors to live.\n", 0, 26);
	liveHigh = getIntInRange("Enter maximum neighbors to live.\n", liveLow, 26);
	maxAge = getIntInRange("Enter maximum age of a point.\n", 1, 10000);
	ageOrInherit = getIntInRange("Enter \"0\" for age coloring, \"1\" for inheritance coloring.\n", 0, 1);

	while(true)
	{
		cout << "Enter path to bmp file (\"0\" for random generation).\n";
		cin >> imagePath;

		if(!cin.good())
		{
			cin.clear();
			cin.ignore(numeric_limits<streamsize>::max(),'\n');
			cout << "Error - enter again.\n";
			continue;
		}
		if(strcmp(imagePath, "0")) // path entered
		{
			if(reader->CanReadFile(imagePath))
			{
				reader->SetFileName(imagePath);
				reader->Update();
				hsv = vtkImageRGBToHSV::New();
				hsv->SetInput(reader->GetOutput());
				hsv->Update();
	
				hsv->GetOutput()->GetDimensions(dim);

				dimX = dim[0];
				dimY = dim[1];

				if(dimX < dimY)
					dimZ = dimX;
				else
					dimZ = dimY;

				useImage = 1;
				break;
			}
			else
			{
				cout << "Cannot open input file - enter again.\n";
				continue;
			}
		}
		else // random generation selected
		{
			//dimX = getIntInRange( "Enter X dimension.\n", 1, 500);
			//dimY = getIntInRange( "Enter Y dimension.\n", 1, 500);
			//dimZ = getIntInRange( "Enter Z dimension.\n", 1, 500);
			useImage = 0;
			break;
		}
	}

	// dimension and timer user inputs are temporarily commented out for our presentation 

	totalSize = dimX*dimY*dimZ;
	layerSize = dimY*dimZ;

	//timer = getIntInRange( "Enter milliseconds between iterations.\n", 10, 10000 );
}

void stepInherit()
{
	// static variables to avoid allocation on every call
	// (just a theory that this improves performance)
	static double neighborColors;
	static int i, j, k, i2, j2, k2, neighborCount, gridOffset,
		offX0, offX1, offY0, offY1, offZ0, offZ1;

	gridOffset = 0;

	for(i = 0; i < dimX; i++)
	{
		// start and stop offsets for neighbor search on x axis
		offX0=-1; offX1=1;
		if(i==0) offX0=0;
		if(i==dimX-1) offX1=0;

		for(j = 0; j < dimY; j++)
		{
			// start and stop offsets for neighbor search on y axis
			offY0=-1; offY1=1;
			if(j==0) offY0=0;
			if(j==dimY-1) offY1=0;

			for(k = 0; k < dimZ; k++)
			{
				// start and stop offsets for neighbor search on z axis
				offZ0=-1; offZ1=1;
				if(k==0) offZ0=0;
				if(k==dimZ-1) offZ1=0;

				neighborCount = 0;
				neighborColors = 0;

				// count living neighbors
				for(i2 = offX0; i2<=offX1; i2++)
					for(j2 = offY0; j2<=offY1; j2++)
						for(k2 = offZ0; k2<=offZ1; k2++)
							if(universe[n][i+i2][j+j2][k+k2])
							{
								neighborCount++;
								neighborColors = neighborColors + scalars[n]->GetValue((i + i2)*layerSize + (j + j2)*dimZ + k + k2);
							}

				// apply game rules here
				if(!universe[n][i][j][k]) // point is dead, check if it will come alive
				{
					if(neighborCount >= bornLow && neighborCount <= bornHigh)
					{
						universe[n2][i][j][k] = 1; // point comes alive
						// color of new point is average of living neighbor colors
						scalars[n2]->SetValue(gridOffset, neighborColors/neighborCount);
					}
					else // point stays dead
					{
						universe[n2][i][j][k] = 0;
						scalars[n2]->SetValue(gridOffset, 1);
					}
				}
				else // point is alive, check if it stays alive or dies
				{
					// decrement count because the point counted itself
					// as a live neighbor in the loop above
					neighborCount--;

					if(neighborCount >= liveLow &&
						neighborCount <= liveHigh &&
						universe[n][i][j][k] <= maxAge)
					{
						universe[n2][i][j][k] = universe[n][i][j][k] + 1; // increment age
						scalars[n2]->SetValue(gridOffset, scalars[n]->GetValue(gridOffset));
					}
					else // point dies
					{
						universe[n2][i][j][k] = 0;
						scalars[n2]->SetValue(gridOffset, 1);
					}
				}
				gridOffset++;
			}
		}
	}
	idata->GetPointData()->SetScalars(scalars[n2]);
}

void stepAge()
{
	// static variables to avoid allocation on every call
	// (just a theory that this improves performance)
	static int i, j, k, i2, j2, k2, neighborCount, gridOffset,
		offX0, offX1, offY0, offY1, offZ0, offZ1;

	gridOffset = 0;

	for(i = 0; i < dimX; i++)
	{
		// start and stop offsets for neighbor search on x axis
		offX0=-1; offX1=1;
		if(i==0) offX0=0;
		if(i==dimX-1) offX1=0;

		for(j = 0; j < dimY; j++)
		{
			// start and stop offsets for neighbor search on y axis
			offY0=-1; offY1=1;
			if(j==0) offY0=0;
			if(j==dimY-1) offY1=0;

			for(k = 0; k < dimZ; k++)
			{
				// start and stop offsets for neighbor search on z axis
				offZ0=-1; offZ1=1;
				if(k==0) offZ0=0;
				if(k==dimZ-1) offZ1=0;

				neighborCount = 0;

				// count living neighbors
				for(i2 = offX0; i2<=offX1; i2++)
					for(j2 = offY0; j2<=offY1; j2++)
						for(k2 = offZ0; k2<=offZ1; k2++)
							if(universe[n][i+i2][j+j2][k+k2])
								neighborCount++;

				// apply game rules here
				if(!universe[n][i][j][k]) // point is dead, check if it will come alive
				{
					if(neighborCount >= bornLow && neighborCount <= bornHigh)
					{
						universe[n2][i][j][k] = 1; // point comes alive
						scalars[0]->SetValue(gridOffset, 0); // scalar 0 -> red
					}
					else // point stays dead
						universe[n2][i][j][k] = 0;
				}
				else // point is alive, check if it stays alive or dies
				{
					// decrement count because the point counted itself
					// as a live neighbor in the loop above
					neighborCount--;

					if(neighborCount >= liveLow &&
						neighborCount <= liveHigh &&
						universe[n][i][j][k] <= maxAge)
					{
						// increment age and set the color scale one fraction higher
						universe[n2][i][j][k] = universe[n][i][j][k] + 1;
						scalars[0]->SetValue(gridOffset, double(universe[n][i][j][k] - 1)/maxAge);
					}
					else // point dies
					{
						universe[n2][i][j][k] = 0;
						scalars[0]->SetValue(gridOffset, 1);
					}
				}
				gridOffset++;
			}
		}
	}
	idata->GetPointData()->SetScalars(scalars[0]);
}

void step()
{
	// switch read/write universe every timestep
	if(n){n = 0; n2 = 1;}
	else{n = 1; n2 = 0;}

	if(ageOrInherit)
		stepInherit();
	else
		stepAge();
}


class vtkSimpleTimerCallback : public vtkCommand
{
	public:
		static vtkSimpleTimerCallback *New()
		{
			vtkSimpleTimerCallback *cb = new vtkSimpleTimerCallback;
			return cb;
		}

		void setTimerToWatch(int timerId)
		{
			timerToWatch = timerId;
		}

		virtual void Execute(vtkObject *caller, unsigned long eventId, void *callData)
		{
			if(vtkCommand::TimerEvent == eventId)
			{
				int timerId = *static_cast<int*>(callData);
				if(timerId != timerToWatch)
					return;
				
				iren->Render();
				step();
				tpoints->Modified();
			}
		}
	private:
		int timerToWatch;
};

int main()
{	
	int i, j, k;
	srand(time(NULL));

	getRules();

	//allocate ptr to [2][dimX][dimY][dimZ] universe
	universe = new int ***[2];
	for(i = 0; i < 2; i++)
	{
		*(universe + i) = new int **[dimX];
		for(j = 0; j < dimX; j++)
		{
			*(*(universe + i) + j) = new int *[dimY];
			for(k = 0; k < dimY; k++)
			{
				*(*(*(universe + i) + j) + k) = new int[dimZ];
			}
		}
	}

	scalars[0] = vtkDoubleArray::New();
	scalars[1] = vtkDoubleArray::New();
	scalars[0]->Allocate(totalSize);
	scalars[1]->Allocate(totalSize);

	if(useImage)
		inputImage();
	else
		randomize();

	idata = vtkImageData::New();
	idata->SetDimensions(dimX, dimY, dimZ);
	idata->SetOrigin(-dimX/2, -dimY/2, -dimZ/2);
	idata->GetPointData()->SetScalars(scalars[0]);

	tpoints = vtkThresholdPoints::New();
	tpoints->SetInput(idata);
	tpoints->ThresholdByLower(RAND_MAX/double(RAND_MAX + 1)); // get living points (scalar < 1)

	idataMapper = vtkPolyDataMapper::New();
	idataMapper->GlobalImmediateModeRenderingOn(); // optimization
	idataMapper->SetInputConnection(tpoints->GetOutputPort());

	idataActor = vtkActor::New();
    idataActor->SetMapper(idataMapper);
	idataActor->GetProperty()->SetPointSize(3);

	renderer = vtkRenderer::New();
	renderer->AddActor(idataActor);
	renderer->SetBackground(0,0,0);
	renderer->GetActiveCamera()->SetPosition(dimX*1.5, dimY*1.5, dimZ*1.5);

	renWin = vtkRenderWindow::New();
    renWin->AddRenderer(renderer);
	renWin->SetSize(600,600);

	iren = vtkRenderWindowInteractor::New();
    iren->SetRenderWindow(renWin);
	iren->Initialize();

	vtkSimpleTimerCallback *cb = vtkSimpleTimerCallback::New();
	iren->AddObserver(vtkCommand::TimerEvent, cb);
	cb->setTimerToWatch(iren->CreateRepeatingTimer(timer));

	iren->Start();

	// clean up

	if(useImage && ageOrInherit)
		lut->Delete();
	reader->Delete();
	if(useImage)
		hsv->Delete();

	scalars[0]->Delete();
	scalars[1]->Delete();
	idata->Delete();
	tpoints->Delete();
	idataMapper->Delete();
	idataActor->Delete();
	renderer->Delete();
	renWin->Delete();
	iren->Delete();

	for(i = 0; i < 2; i++)
	{
		for(j = 0; j < dimX; j++)
		{
			for(k = 0; k < dimY; k++)
			{
				delete[] *(*(*(universe + i) + j) + k);
			}
			delete[] *(*(universe + i) + j);
		}
		delete[] *(universe + i);
	}
	delete[] universe;

	return 0;
}